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Abstract. - We compare the vibrational properties of model Si02 glasses generated by 
molecular-dynamics simulations using the effective force field of van Beest et al. (BKS) with 
those obtained when the BKS structure is relaxed using an ab initio calculation in the frame- 
work of the density functional theory. We find that this relaxation significantly improves the 
agreement of the density of states with the experimental result. For frequencies between 14 
and 26 THz the nature of the vibrational modes as determined from the BKS model is very 
different from the one from the ab initio calculation, showing that the interpretation of the 
vibrational spectra in terms of calculations using effective potentials can be very misleading. 

Motivation. Understanding the microscopic properties of vibrational excitations in 
disordered systems is a long standing challenge in basic physics as well as in material science 
since the lack of positional order makes both the experimental study and the theoretical 
interpretations of the results very difficult. For instance the mechanism leading to the existence 
the so-called boson peak present in many glasses is the subject of a long standing debate, and 
the reason for the presence of the Di and D2 lines in the Raman spectra of amorphous SiC>2 
has remained unclear for a long time . 

In principle molecular dynamics (MD) computer simulations overcome these difficulties 
since one has direct access to all the necessary microscopic information. Therefore in recent 
years many studies of this kind have been carried out with the aim of shedding some light on 
the nature of these vibrational excitations, in particular for the case of silica, the paradigm 
of network forming glasses [p|-jl9||. Due to the large computational costs of such simulations 
the vast majority of them were done with effective classical force fields, i.e. potentials which 
were optimized to reproduce certain (somewhat arbitrarily chosen) experimental features of 
SiC>2 • It is clear that the reliability of the results of these investigations depends crucially on 
whether or not the interactions used are sufficiently accurate to allow a faithful description 
of the real material. Hence a considerable effort has been made to check that the models 
used do reproduce the salient structural and dynamical features of real silica. In these studies 
it has been shown that the classical force fields employed are indeed able to give a good 
description of quantities like the structure factor, the diffusion constant or viscosity, etc. jTJJ 
pfj|-p^] , so their use in investigations of the vibrational properties also appears a reasonable 
undertaking. Nevertheless, it was found that certain features of the vibrational density of 
states (DOS) are not well reproduced by these models and more sophisticated calculations, 
such as the numerical studies based on first-principles, seem to be required p5fl , in agreement 
with conclusions drawn for the case of liquids [ p6[ . Pasquarello et al. showed that using an ab 
initio scheme it is possible to reproduce many structural, electronic and vibrational properties 
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of real silica ||.[U]. [12|,[1J|. However this approach suffers from its heavy computational cost 
which restricts this type of calculation to the study of very small systems with a rather low 
statistical accuracy. 

One possibility to overcome this limitation, at least partially, is to use a combined approach 
which consists in generating a glass using an effective potential and subsequently refining the 
structure obtained by means of first-principles j27j. In previous work we showed that the 
structure of vitreous SiC>2 generated using the effective force field by van Beest, Kramer and 
van Santen (BKS) |28| is only modified weakly by a first-principles calculation, thus validating 
the structural model generated with this potential. 

In contrast to this we show in this letter that the DOS of a Si02 glass generated by classical 
MD simulations using the BKS potential is strongly modified by using an ab initio treatment 
of the forces, and that this treatment leads to a much better agreement with experimental 
results. Moreover, in a large frequency range, the nature of the excitations as determined from 
the effective potential differs significantly from the one determined from the ab initio forces 
thus raising doubts as to the detailed analysis of the nature of the vibrational excitations 
determined from the BKS force field. 

Simulation details. - Molecular-dynamics simulations were done using the BKS potential 
on systems containing 26 Si02 units at the experimental density (2.2 g/cm 3 ). For this we 
used the velocity form of the Verlet algorithm with a time step of f.63 fs. Three different 
samples were generated by quenching liquids well-equilibrated at 3500 K to 300 K, using three 
different cooling rates: 5-10 12 K/s, 3 TO 11 K/s, and 7-10 10 K/s. The glasses obtained this way 
(which are non-equilibrium structures) were annealing for 70 ps at 300 K, and subsequently 
quenched to K, at which their dynamical matrices were evaluated and diagonalized in order 
to obtain the vibrational frequencies and the corresponding (normalized) eigenmodes. In 
parallel the final atomic coordinates and velocities after the annealing at 300 K were used as 
initial coordinates and velocities for short (« 0.12 ps) ab initio molecular-dynamics simulations 
of the Car-Parrinello type p9[ , using the CPMD code J3(|. The technical details of these 
simulations were identical to the ones described in Ref. |27[ . At the end of these simulations 
the structures of the three glasses were relaxed to K and the dynamical matrices were 
computed by evaluating the second derivatives of the total energy with respect to atomic 
displacements by taking finite differences of the atomic forces. Subsequently the vibrational 
frequencies and the corresponding eigenmodes were obtained from these matrices. Hence we 
obtained g(co), the true DOS for this system. Note that although the cooling rates are high 
and the system size is small, the DOS depends only weakly on these parameters |pof . 

In the following, the quantities computed by means of classical molecular-dynamics sim- 
ulations using the BKS potential and the CPMD code will be labeled "BKS" and "CP", 
respectively. 

Results. - Since we found that the DOS from the three different cooling rates are identical 
to within the statistical error - which is relatively large due to the small system size -, we 
decided to treat the three glasses as independent statistical samples and analyzed the three 
sets of vibrational frequencies/modes together. The resulting vibrational DOS was used to 
compute an effective neutron scattering cross section G(lu) ~ C{uj)g(uj). This was done by 
using the incoherent approximation and by calculating C(oj) as in Ref. Jlgj . We note that the 
correction functions C(lo) for the BKS and CP are very similar and hence differences in the 
respective G(tu) are mainly due to differences in the respective g{u>). In Fig. [I] we compare the 
G(w) obtained and we see that at intermediate frequencies the two curves are very different. 
In particular we see that the CP curve has a pronounced peak at around 12 THz and a 
smaller one at around 24 THz. Finally there is a small peak at 18 THz, the so-called D2 



Magali Benoit and Walter Kob: 



3 




Fig. 1 - Frequency dependence of the effective density of states G(ui) calculated from ab initio (solid 
line) and classical (dashed line) molecular dynamics simulations and compared to neutron scattering 
experiments from Ref. j^l] (filled circles). The calculated G(cj) have been obtained using the exper- 
imental values for the neutron scattering length factors bsi = 4.149 fm and 6o = 5.803 fm and a 
Gaussian broadening of width 2a = 1.05 THz. 



line, which is due to a ring of size three. Overall the CP curve is in very good agreement 
with previous investigations Q . All these features are missing in the BKS curve, despite the 
good agreement between CPMD and BKS with regard to structural properties. Also included 
is the result of a neutron scattering experiment by Carpenter and Price pl[ and from the 
reasonable agreement between this curve and the one from the CP calculation we conclude that 
the latter is reliable. [ Note that i) there is no fit parameter whatsoever, and ii) the lack of a 
small peak at around 4 THz in the experimental data is related to the insufficient experimental 
resolution ]. Hence we conclude from this figure that the DOS as calculated from the BKS 



model is not very reliable at intermediate frequencies, in agreement with Ref. 11 1. Note that 
similar discrepancies between experiments and simulations with various effective interactions 
have already been observed in previous studies [pl],[il||i^,|2rj|,|2^| and hence we conclude that 
many other types of force fields also lead to a density of states which is not trustworthy and 
that most probably the conclusions drawn in this paper hold for these other potentials as well. 
However, the good agreement between the CP and the experimental DOS clearly shows that 
a simple refinement of the BKS model glass by a first-principles calculation is sufficient to 
significantly improve the vibrational density of states. 

Having shown that the effective DOS from the BKS differs significantly from the one from 
the CP we now investigate this difference in more detail by comparing the corresponding 
eigenmodes. Since the number of eigenmodes for the BKS and CP cases is the same, Fig. [j] 
implies that there is a considerable reshuffling of the modes from one frequency to another 
if the force field is switched. In order to check whether at a given frequency 
1, . . . , 3 x 3A^, and N is the number of atoms in each of the three samples) the nature of a 
given mode e BKS (w BKS ) from the BKS system is similar to one of the CP modes we calculated 
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Fig. 2 - Contour plot of ,u£ F ), 

the projection of the BKS eigenmodes onto the CP modes 
(see Eq. ([!])). The solid lines are the DOS g{ui) as determined from the BKS and the CP glasses 
respectively, scaled in order to fit into the graph. [Contour levels: 0.06, 0.1, 0.15, 0.2, 0.3, 0.5, 0.8 
1.0] 



the projections of the latter onto the former modes: 

i/K BKS ,0 = |e BKS (^ BKS )-e CP «)| with „ >M =l,...,3xaiV. (1) 

Here e CP (w^ p ) is the eigenvector for the CP system at frequency w^ p . With this definition, 
H varies between and 1, since |e BKS (u; BKS )| = |e CP (u;^ p )| = 1. If the nature of the modes in 
the BKS system were identical to the ones of the CP system, H would be a line of 6— functions 
along the diagonal in the (w BKS ,w^ p ) plane. The frequency dependence of H(iv BKS , u>^ p ) is 
shown in Fig. || as a contour plot. 

This plot shows a well-defined ridge close to the diagonal of the graph up to w 7 THz 
and between rj 28 and 39 THz. This means that the overlap between the BKS and the CP 
vibrational modes is significant in these frequency ranges, i.e. that at a given frequency the 
nature of the mode is very similar. For 7 THz < ui CF < 13 THz we still find a well defined 
ridge, but its location is above the diagonal. Hence we see that although the BKS modes can 
still be well represented by the CP ones, they are shifted to slightly higher frequencies. Since 
the interval 7 THz < lu bks < 21 THz is compressed into the range 7 THz < lu cp < 14 THz 
this leads to a significant larger g CF (tu) in this range, in agreement with Fig. [l]. 

In the range 16 THz < uj bks < 28 THz we find that the BKS modes have a significant 
overlap with CP modes that cover a large range in w CP . This implies that in this frequency 
range the vibrational dynamics of the BKS system is no longer realistic and that hence care 
must be taken if one draws conclusions from the analysis of the modes in this range. 

In order to investigate whether in this frequency range the CP eigenmodes can at least 
be written as a sum of only a few BKS modes, we ordered for each frequency oj^ p the over- 
laps ff(w BK VCP) in descending order, i.e. [H (uj bks , uj^ p )} 2 > [H(tu BKS ,w° F )] 2 > . . ., and 
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Fig. 3 - Circles: Sum of the J largest overlaps between the BKS and CP modes, as defined in Eq. (g) , 
for each CP eigenfrequency wj; p . Bold solid line: The CP-DOS (scaled in order to fit in the graph). 



calculated the sum 

^K CP )-E[^K B ' KS ^, CP )] 2 • (2) 
i/'=i 

(Note that the sum of such eigenmodes is in general no longer an eigenmode.) The resulting 
Sj as a function of the CP vibrational frequencies are presented in Fig. || for J— I, 2, 10, 
and 20. This figure demonstrates that in general the CP modes are not well represented until 
the sum extends over at least 10 BKS eigenmodes. The only exceptions are modes at very 
low frequency (< 3 THz) and around 27, 29, and 38 THz. These latter modes show already 
a prominent peak in Si(u>^ p ) and are at band edges and hence are localized. Note that even 
for 3 THz < u CP < 13 THz, i.e. the frequency range where we find a strong correlation in 
Fig. ||, the BKS modes do not describe well the CP modes. In the range 13 THz < lu cp < 26 
THz one needs on the order of 20 BKS modes in order to describe a CP mode, which shows 
that the former have not much in common with the latter. 

It is well known that the vibrational excitations of silica include modes that are very 
localized as well as ones that are collective. In order to see whether there is a difference in 
the ability of the BKS potential to describe the one or the other type we have defined the 
following quantities: 

X7(^ P ) =EE |e^K CP )f and X*™(u>™ s ) = £ £ |e^ s (^ KS )| 2 (3) 

ct' — l i—1 a' — l 2—1 

where a' is the atomic index, i = x,y,z and the prime in the sum (a') indicates that we have 
made the ordering: J2^=i \ e i,i\ 2 > Si=i l e 2,i| 2 > • ' • Thus X^ p (uj^ p ) is the sum of the J 
largest atomic displacements for that eigenvector. The resulting X BKS and X CP are shown 
in Figure ||, for values of J=l, 4, 20 and 50. (Note that for the sake of clarity the quantity is 
only shown for one sample. The other two samples look qualitatively similar.) 
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Fig. 4 - Sum of the J largest components of an eigenvector with frequency ui (see Eq. (|3|)) for J = 
1, 4, 20 and 50. Left panel: X^; Right panel: X^ KS . 



Using this quantity, we can estimate the number of atoms participating in a given vibra- 
tional mode. Firstly, we note that most of the modes with frequencies below ps 22 THz have a 
"collective" nature since around 20 eigenvectors are needed to make 60 % of X, independent 
of the model used. Exceptions are CP modes at 14 and 18 THz which show a pronounced 
peak with a height larger than 0.4 even for J = 4. These modes correspond to the Di and D2 
lines in the Raman spectrum of vitreous Si02 that have been attributed to breathing modes 
of 3 and 4-membered rings jll|[l!|. From the curves it becomes clear that these modes are not 
present in the BKS description, thus raising the problem of the interpretation of the Raman 
spectrum of glassy SiC-2 using this model potential (see also p5|). For frequencies higher 
than « 22 THz the nature of the modes becomes more localized for the case of the CP model 
as well as the BKS model. However, we note a clear difference in the Xj KS and Xj p between 
22 and 29 THz, where the BKS modes seem to be more collective than the CP ones. Finally 
we mention that one sees a strong localization of the modes near the gap of the DOS and at 
around 39 THz, in agreement with the comment made in context of Fig. ||. 

Conclusion. We have shown that for the case of silica the DOS as predicted by an 
effective force field can be improved considerably by relaxing the configuration by means of an 
ab initio simulation, without a significant change in the structural properties. For intermediate 
frequencies the nature of the modes in the BKS system is very different from the one in the 
CP system. In particular the effective force field is not able to reproduce the nature of the 
modes attributed to the Di and D2 lines in the Raman spectra. In view of the fact that 
the BKS model was obtained by an ab initio calculation of a single tetrahedron and a lattice 
dynamic simulation of a crystal to optimize the elastic constants psj ] , it is not surprising that 
the model does well at very high and very low frequencies, but is not reliable at intermediate 
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frequencies. Hence we conclude that despite the reliability of effective potentials with regard 
to structural properties, it might be that if one wants to investigate vibrational features in 
the whole frequency range one would need an effective potential which is more accurate, or a 
full ab initio calculation. 

* * * 

We thank J. Horbach, S. Ispas and R. Jullien for useful discussions. Computations have 
been performed on the IBM SP3 of the national computer center CINES in Montpcllicr, 
France. 
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